home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / csim / source.lha / source / C++SIM / Random.cc < prev    next >
C/C++ Source or Header  |  1993-06-14  |  4KB  |  199 lines
  1. /*
  2.  * Copyright (C) 1993
  3.  *
  4.  * Department of Computing Science,
  5.  * The University,
  6.  * Newcastle upon Tyne,
  7.  * UK.
  8.  */
  9.  
  10. #include <iostream.h>
  11. #include <math.h>
  12.  
  13. #ifndef RANDOM_H_
  14. #include "Random.h"
  15. #endif
  16.  
  17.  
  18. RandomStream::RandomStream (long MGSeed, long LCGSeed)
  19. {
  20.     // Clean up input parameters
  21.     if ((MGSeed&1) == 0) MGSeed--;
  22.     if (MGSeed<0) MGSeed = -MGSeed;
  23.     if (LCGSeed<0) LCGSeed = -LCGSeed;
  24.  
  25.     // Initialise state
  26.     MSeed = MGSeed;
  27.     LSeed = LCGSeed;
  28.  
  29.     for (int i=0; i< (sizeof(series)/sizeof(double)); i++)
  30.     series[i] = MGen();
  31. }
  32.  
  33. double RandomStream::MGen ()
  34. {
  35.     // A multiplicative generator, courtesy I. Mitrani 1992,
  36.     // private correspondence
  37.     // Y[i+1] = Y[i] * 5^5 mod 2^26
  38.     // period is 2^24, initial seed must be odd
  39.  
  40.     const long int two2the26th = 67108864;    // 2**26
  41.  
  42.     MSeed = (MSeed * 25) % two2the26th;
  43.     MSeed = (MSeed * 25) % two2the26th;
  44.     MSeed = (MSeed * 5) % two2the26th;
  45.  
  46.     return (double) MSeed / (double) two2the26th;
  47. }
  48.  
  49. double RandomStream::Uniform () 
  50. {
  51.     // A linear congruential generator based on the algorithm from
  52.     // "Algorithms", R. Sedgewick, Addison-Wesley, Reading MA, 1983.
  53.     // pp. 36-38.
  54.     const long m=100000000;
  55.     const long b=31415821;
  56.     const long m1=10000;
  57.  
  58.     // Do the multiplication in pieces to avoid overflow
  59.     long    p0 = LSeed%m1,
  60.         p1 = LSeed/m1,
  61.         q0 = b%m1,
  62.         q1 = b/m1;
  63.  
  64.     LSeed = (((((p0*q1+p1*q0)%m1)*m1+p0*q0)%m) + 1) % m;
  65.  
  66.     // The results of the LC generator are shuffled with
  67.     // the multiplicative generator as suggested by
  68.     // Maclaren and Marsaglia (See Knuth Vol2, Seminumerical Algorithms)
  69.  
  70.     long choose = LSeed % (sizeof(series)/sizeof(double));
  71.  
  72.     double result = series[choose];
  73.     series[choose] =  MGen();
  74.  
  75.     return result;
  76. }
  77.  
  78. double RandomStream::Error ()
  79. {
  80.     const long r=100;
  81.     const long N=100*r;
  82.     long i, f[r];
  83.     for (i=0; i<r; i++) f[i]=0;
  84.     for (i=0; i<N; i++) f[(int) (Uniform()*r)]++;
  85.     long t=0;
  86.     for (i=0; i<r; i++) t += f[i]*f[i];
  87.     double rt = (double) r*t;
  88.     double rtN = rt / (double) N - (double) N;
  89.     return 1.0 - (rtN / r);
  90. }
  91.  
  92. UniformStream::UniformStream (double l, double h, int StreamSelect)
  93. {
  94.     lo = l;
  95.     hi = h;
  96.     range = hi-lo;
  97.     for (int ss=0; ss<StreamSelect*1000; ss++) (void) Uniform();
  98. }
  99.  
  100. double UniformStream::operator() ()
  101. {
  102.     return lo+(range*Uniform());
  103. }
  104.     
  105. Draw::Draw(double p, int StreamSelect) : s(0,1)
  106. {
  107.     prob = p;
  108.     for (int ss=0; ss<StreamSelect*1000; ss++) (void) s();
  109. }
  110.  
  111. boolean Draw::operator() ()
  112. {
  113.     double x = s();
  114.  
  115.     if (x >= prob)
  116.     return true;
  117.     else
  118.     return false;
  119. }
  120.  
  121. ExponentialStream::ExponentialStream (double m, int StreamSelect)
  122. {
  123.     Mean = m;
  124.     for (int ss=0; ss<StreamSelect*1000; ss++) (void) Uniform();
  125. }
  126.  
  127. double ExponentialStream::operator() ()
  128. {
  129.     return -Mean*log(Uniform());
  130. }
  131.  
  132. ErlangStream::ErlangStream (double m, double s, int StreamSelect)
  133. {
  134.     Mean = m;
  135.     StandardDeviation = s;
  136.  
  137.     double z = Mean/StandardDeviation;
  138.     k = (long) (z*z);
  139.     for (int ss=0; ss<StreamSelect*1000; ss++) (void) Uniform();
  140. }
  141.  
  142. double ErlangStream::operator() ()
  143. {
  144.     double z=1.0;
  145.     for (int i=0; i<k; i++) z*=Uniform();
  146.     return -(Mean/k)*log(z);
  147. }
  148.  
  149. HyperExponentialStream::HyperExponentialStream (double m, double s, int StreamSelect)
  150. {
  151.     Mean = m;
  152.     StandardDeviation = s;
  153.     double cv,z;
  154.     cv=StandardDeviation/Mean;
  155.     z = cv*cv;
  156.     p = 0.5*(1.0-sqrt((z-1.0)/(z+1.0)));
  157.     for (int ss=0; ss<StreamSelect*1000; ss++) (void) Uniform();
  158. }
  159.  
  160. double HyperExponentialStream::operator() ()
  161. {
  162.     double z = (Uniform()>p) ? Mean/(1.0-p) : Mean/p;
  163.     return -0.5*z*log(Uniform());
  164. }
  165.  
  166. NormalStream::NormalStream (double m, double s, int StreamSelect)
  167. {
  168.     Mean = m;
  169.     StandardDeviation = s;
  170.     z = 0.0;
  171.     for (int ss=0; ss<StreamSelect*1000; ss++) (void) Uniform();
  172. }
  173.  
  174. double NormalStream::operator() ()
  175. {
  176.     // Use the polar method, due to Box, Muller and Marsaglia
  177.     // Taken from Seminumerical Algorithms, Knuth, Addison-Wesley, p.117
  178.  
  179.     double X2;
  180.  
  181.     if (z!=0.0) {
  182.     X2 = z;
  183.     z = 0.0;
  184.     } else {
  185.     double S, v1, v2;
  186.     do {
  187.         v1 = 2.0*Uniform()-1.0;
  188.         v2 = 2.0*Uniform()-1.0;
  189.         S = v1*v1 + v2*v2;
  190.     } while (S>=1.0);
  191.  
  192.     S = sqrt((-2.0*log(S))/S);
  193.     X2 = v1*S;
  194.     z  = v2*S;
  195.     }
  196.  
  197.     return Mean + X2*StandardDeviation;
  198. }
  199.